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Abstract 

An algorithm for calculating the spectral intensity of radiation due to the co- 
herent addition of many particles with arbitrary trajectories is described. Di- 
rect numerical integration of the Lienard-Wiechert potentials, in the far-field, 
for extremely high photon energies and many particles is made computation- 
ally feasible by a mixed analytic and numerical method. Exact integrals of 
spectral intensity are made between discretely sampled trajectories, by as- 
suming the space-time four-vector is a quadratic function of proper time. The 
integral Fourier transform of the trajectory with respect to time, the modu- 
lus squared of which comprises the spectral intensity, can then be formed by 
piecewise summation of exact integrals between discrete points. Because of 
this, the calculation is not restricted by discrete sampling bandwidth theory, 
and hence for smooth trajectories, time-steps many orders larger than the 
inverse of the frequency of interest can be taken. 
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1. Introduction 

Radiation from synchrotrons is a well developed field, and a number of 
numerical methods for calculating the radiation spectrum exist [H El |3l HI [5] . 
Laser driven sources of radiation have also triggered interest in measurements 
illZllHllSlIinillllllSllISlEj and numerical calculations [ISl HSl CZl HE] of 
radiation, particularly in the area of laser-plasma particle acceleration. As 
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sources of high-energy particle beams and radiation, laser-plasma based tech- 
niques may be used for a large range of future applications. Ultrafast X-ray 
sources would be useful in, for example, time resolved diffraction, medical 
imaging, spectroscopy and microscopy of transient physical, chemical, or bio- 
logical phenomena. Laser wakefield acceleration [12] of high energy electrons 
beams has recently successfully demonstrated the production of GeV peak 
energy electron beams |20j , and has become a highly cited field of research. 
Such beams could be used directly for radiotherapy or radiographic imaging, 
or alternatively can be converted into fs duration, high brightness sources 
of x-rays. Numerical calculation of the x-ray spectrum is therefore of in- 
terest. Scattering of laser pulses from relativistic electron beams is another 
area which may require well characterized angularly resolved spectra from a 
realistic bunch, for comparison with experiment. 

In general, the trajectories of particles in laser driven experiments are 
quite complicated; wakefields, electron beam interactions with intense laser 
fields, electron orbits in laser generated channels and in laser-solid interac- 
tions all represent sources of radiation. The radiation fields can be explicitly 
calculated by a fast Fourier transform method or by finite differencing of 
the Lienard-Wiechert fields in time, but these are computationally intensive 
processes due to the constraints of the Whittaker-Shannon-Nyquist sampling 
theorem [211 [22] • This essentially states that for discretely sampled data, fre- 
quencies higher than half the sampling frequency are aliased to frequencies 
lower than half the sampling frequency. Since the spectral power as a func- 
tion of frequency emitted is essentially equivalent to a Fourier transform, the 
highest resolvable frequency in the spectrum is constrained. However, for 
relativistic particles radiation can be produced at much higher frequencies 
than the actual time-scale corresponding to the change in momentum of the 
particles, because the radiation co-propagates with the particle. 

Here, an algorithm is developed which uses a combination of numerical 
and analytic methods to integrate the Lienard-Wiechert fields, for arbitrary 
trajectories of many charged particles, to frequencies greatly exceeding the 
Nyquist frequency, I'm- This means that the sampling rate to accurately re- 
produce a particular spectrum can be many orders of magnitude lower than 
I/vni and therefore faster to solve numerically. This is particularly rele- 
vant to calculations of radiation from betatron oscillations in laser wakefield 
accelerators and beam-laser interactions, but the technique is applicable to 
numerous other areas of physics. 
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2. The numerical method 



The spectral intensity of radiation emitted by a number Np of accelerat- 
ing point charges, with the jth particle at position and with normalized 
velocity Pj = v^/c, can be expressed, in the far-field, as: 
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(1) 



where the unit vector § is in the direction of observation, at a distance far 
compared with the scale of the emission region. This can be written alter- 
natively in terms of proper time, r: 
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where the four-wave-vector = uj{^, s/c}, and Vj is the momentum part of 
the jth particle's four-velocity defined as: 



dXj 



(3) 



One of the advantages of using proper time rather than 'laboratory' time for 
numerical calculations is that for a uniform step finite differencing scheme, 
the time resolution is effectively adaptive; as the particle gains inertia and 
is therefore accelerated at a decreased rate for a similar force, the laboratory 
time step-size increases. To numerically integrate the equations of motions for 
charged particles, both and v have to be recorded at a number of discrete 
points. To then perform the spectral integration numerically, a 'zeroth order' 
model would be to reduce equation [2] to a summation over finite differenced 
points: 
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EE 

j=l n=0 

The problem with this method is that frequencies higher than half the sam- 
pling rate are aliased to lower frequencies, and therefore an upper limit is 
put on the maximum frequency that can be effectively resolved to ~ 1/Ar 
[22] • For attempts to simulate high energy photons from laser interac- 
tions, this can be computationally prohibitive. However, the motion of the 
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particles which lead to such high-energy photons generally consists of changes 
on timescales much larger than the radiation frequencies produced. Hence, a 
different approach is taken, which is to assume that the motion of a particle 
between time-steps can be approximated by an interpolating function. The 
spectral integral can be expressed as a summation over analytic integrals 
between each time-step: 

5^5^sx / v,(r)e--."«dr . (5) 



(Pi jUqC^C 2 



=1 n=0 



Here the notation used is that Greek character sub/super-scripts represents 
components of 4- vector quantities, subscript Roman characters denote either 
particle number, j, or (proper) time-step, n, and bold font represents 3- vector 
quantities. Proper time is discretized into steps of size Ar. An exact analytic 
solution can then be employed to calculate the sub-integral between time- 
steps, using an interpolating function for Vj and between discrete time 
steps. A quadratic interpolation for can accurately model both linear and 
harmonic accelerations, and a linear interpolation for the velocity four-vector 
is consistent. Thus, the four- velocity and four-displacement at a time r — r„ 
are approximated by: 

X (r) = a;Q„ + X]^„(r — r„) + a;2„(T — r„) , 

= Voj,„ + vij-„(r - r„) , (6) 

where Xq„ = x°'(r„), Voj,n = Vj(r„), and Vij„, x^^ and x^n are interpolation 
coefficients. Each integral over time in equation [5] can then be broken up 
into a series of analytic integrals between time steps from r = riAr to r = 
{n + l)Ar: 
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dujdVt IGtt^ 
where, with quadratic interpolation: 

Jr,i-AT/2 

and the space and velocity fourvectors are assumed real so that X*„(co') = 
Xj^„(— cj). A change of variables in the integral leads to: 

J-At/2 
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If each particle trajectory is accurately described by the interpolation func- 
tion between grid-points, then the total integral is exactly solved. This 
is the key difference which allows calculation of the spectrum to far be- 
yond the Nyquist frequency corresponding to the time step between known 
position and moment values. The integral with respect to r, Xj,j(a;) = 
exp{—iKaXQj,^)Ij^n{^^), cau be split into real and imaginary parts: 



/n.T/ z 
(V0j,n + Vij- „r) COS [KaXlj „T + K^x'^j^rj'^) dr , (10) 
■Ar/2 

/Ar/2 
(voj> + vij-„r) sin [k^xI^^t + k,^x"^^t'^) dr . (11) 
-Ar/2 

It can be shown that these integrals have the solutions: 

r *+[c(0+)-c(0_)] 

5J(X,,„) = <^ + [^(6+) - 5(0_)] (12) 

^A2i,n y ^ 2vij-„ [sin$+ - sin<l)_] 

r *+[5(e+)-5(0_)] 

S>(J,,„) = ±- I - [C(0+) - C(0_)] (13) 

^X2j,n y _ 2Vij-„[cOS$+ -COS$_] 

where C{x) = f^^ cos^nt"^ /2)dt and S{x) = sm^irt'^ / 2) dt are the Fresnel 
integrals, and: 

Xlj,n ± X2j,nAr 

'-'i - /o — ' y^^) 

*+ = (2X2i,nVoj> - Xli,nVlj>)cOS r^^^) , (15) 

= J (2X2i,nVoj,„ - Xij,nVij>)sin r^i^^ , (16) 
Ar^ Ar 

^± = -^X2j,n±—Xlj,n, (17) 

Xlj,n ^a-^lj,n i (l^) 

X2j,n = ^a^2j,n • (19) 

Here, the Fresnel integrals are solved numerically by using the power se- 
ries and continued fraction expressions in Numerical Recipes in C++ 
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Although the benefit gained in using less time steps by this method far 
outweighs the cost of calculating these functions, this is a computationally 
expensive process, and if the use of Fresnel integrals can be avoided it would 
be beneficial. If the second exponent in the integral is small, it is appropriate 
to Taylor expand the exponential function and truncate at order r^, hence: 
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The integrals when x'^j^n'^'^ is small are: 
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where sinc(x) is the unnormalized cardinal sine function, sinc(x) = sm{x)/x. 
The Taylor expanded solution for the integral Ij^n is significantly faster to 
solve than the Fresnel integral expressions, so a condition statement can 
switch between solutions depending on the size of X2j,n^T'^- For practical 
purposes, the expansion solution is necessary because the exact solution, 
when numerically solved, diverges for small X2j,n^'^'^ due to floating point 
truncation error, the accuracy of the trigonometric functions and the inverse 
X2j>Ar^ relationship. Some threshold value, T, for X2j,n^T'^ can be used to 
choose between the two models. This must be between 1 and 0, as the Taylor 
expansion is only valid for X2j,n^'T'^ < 1- For values less than the threshold, 
X2j,n^'^'^ < T, the Taylor expansion solution is used and for values greater 
than threshold, X2j.n^T'^ > T, the exact solution is used. In addition, long 
double precision floating point numbers are necessary for sufficiently accurate 
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calculation of the functions. A good threshold value is T = 10~^, as a balance 
between the accuracy of the solution and the speed of the algorithm. 

Figure [T] shows calculations of integral Tj_„ for dr = 1, foj> = 1, vij,n = 2, 
Xij,n = 1 and varying X2j,n- (a) Calculation of integral using Fresnel integral 
form, (b) Calculation using Taylor series expansion form of integral, (c) 
Calculation using finite differenced integral with 10^ steps and Simpson's 
rule. It is clear in this specific example that the Fresnel integral solution is 
correct for all but the smallest values of X2j,n- For these values, since the 



equations 12 and 13 involve the subtraction of large terms and division by 
X2j,n, the truncation errors are amplified to significant values. The Taylor 
expanded solution does not contain a term l/x2j,n, and therefore approaches 
the correct result for X2j,n < 1. 
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Figure 1: Various model calculations of integral X,-.„ for dr — 1, woj.n — 1, vij^n = 2, 
Xlj.n — 1 and varying X2j,n- (a) Calculation of integral using fresnel integral form, (b) 
Calculation using Taylor series expansion form of integral, (c) Calculation using finite 
differenced integral with 10^ steps and Simpson's rule. 



The radiated spectral intensity is then given by: 
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The radiation would normally be considered in a spherical polar coordinate 
system, {r, ^, 0}, where 9 is the azimuth and is the polar angle. For sim- 
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plicity, a cartesian coordinate system, x, ?/, can be chosen such that the 
radiation is calculated at an angle 9 with respect to the z axis in the y — z 
plane. To observe radiation in a particular direction in the coordinate 
system can be rotated about the z axis. Defining x,y, z components of com- 
binations of integrals Xj 

f Np ^ 

r Np Nr -\ 

^iS,,y,.) = X] 5^ (53 (J,,0 cosKx-g + 3? (J,,0 sinK^^^^^^ (2^) 
the radiated spectral intensity is: 



A^ = T^^' +(5^(5.)cos^^-5R(5.)sin^f 
^'"^^ ' +($>(4)cos^^-^>(^.)sin^^f 



(29) 



3. Numerical benchmarks 

3.1. Synchrotron radiation 

A test for the algorithm for relativistic motion is reproducing the well 
known continuum spectrum from a single ultra-relativistic electron rotating 
in a magnetic field without losing energy. Mathematically this means no radi- 
ation damping term in the force equation, physically this can be accelerating 
electric fields that compensate for the energy losses. An analytic expression 
for the spectrum is the sum of modified Bessel functions of the second kind. 
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where ^ = up/Scj^ (1 + Q'^^'^'f^'^ and p is the radius of curvature. Figure [2] 
shows calculated synchrotron spectra (left) and the percentage error relative 
to the analytic solution (right) for different numerical integration techniques. 
The relative percentage error is defined as the difference between the ana- 
lytic and numerical solutions divided by the maximum value of the analytic 
solution. 
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Figure 2: Calculated synchrotron spectra (left) and the percentage error relative to the an- 
alytic solution (right) for different numerical integration techniques. The electron Lorentz 
factor is 7 = 1000, and the magnetic field strength, Bq, is such that ei?o/"^e'^o = 1- The 
time-step is ojqAt = it x 10~'^. (a-b) Calculated spectrum and error for finite differenced 
numerical integration, (c-d) Calculated spectrum and error for summation of exact inte- 
grals using first-order interpolation of position, (e-f) Calculated spectrum and error for 
summation of exact integrals using second-order interpolation of position. In the left-hand 
images, the red dashed line corresponds to the analytic solution, and the blue solid line 
corresponds to the numerically calculated solution. 
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The electron Lorentz factor is 7 = 1000, and the magnetic field strength, 
B, is such that eB/rrieUJo = 1. The time-step is uqAt = tt x 10~^, and the 
maximum frequency calculated is IO'^cuq. Note that this means that the max- 
imum frequency resolved is tt x 10^ times the inverse of the time-step size, 
(a-b) show the calculated spectrum and error for finite differenced numerical 
integration, (c-d) show the spectrum and error for the summation of exact 
integrals using first-order interpolation of position only and (e-f) show the 
spectrum and error for the summation of exact integrals using second-order 
interpolation of position. In the left-hand images, the red dashed line cor- 
responds to the analytic solution, and the blue solid line corresponds to the 
numerically calculated solution. It is clear from (a-b) that using the simple 
finite differencing method (equivalent to a discrete Fourier transform) that 
the resulting spectrum is of no use whatsoever with this particular time-step. 
The error grows with frequency to be 12 orders of magnitude larger than the 
maximum of the spectral intensity that the method is trying to reproduce. 
Decreasing the time-step size eventually yields an error smaller than 100%, 
but to produce accurate spectra, uqAt ^ 1. This is a very limiting factor 
in calculations of this kind. 

By performing a first-order interpolation, as in figure [2] (c-d), which corre- 
sponds to using the real part of the Taylor expanded solution only, equation 



21 the calculated spectrum is clearly now representative of the analytic spec- 



trum, albeit with an error of up to 10%. The simplicity of this algorithm 
(the only non algebraic operation involving the calculation of a single sinusoid 
for each time-step) makes it fast and easy to implement. However, beyond 
lO^cuo, the error in the numerical solution starts to grow. Using second-order 
interpolation, (e-f), yields an error that is at most more than an order of 
magnitude smaller than the maximum error in the first-order interpolation 
calculation. Importantly, it is free of spurious oscillations which could lead 
to misinterpretation of more complex spectra. In addition, the error remains 
less than the maximum error shown in figure up to 5 x 10^a;o, which 
corresponds to a frequency of u > 10^/Ar. Despite being a more complex 
algorithm, the second-order method greatly speeds up a spectral calculation 
due to a larger time-step being allowed. Provided the motion of a particle 
can be well-described by piecewise quadratic functions (i.e. cubic or higher 
terms would be small), the method will produce accurate spectra. 

In figure |3j angularly resolved polarization components of synchrotron 
spectra for the same parameters are shown. In this case, the frequency scale 
is shown as an energy scale corresponding to a magnetic field of Bq = 13 T, 
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Figure 3: Calculated angularly resolved polarization components of synchrotron spectra 
for an electron Lorentz factor of 7 = 1000, and a magnetic field strength of Bq = 13 T. 
The time-step is uiqAt — tt x 10^"*. (left) uj_ is the radiated spectral power polarized in 
the plane of the motion of the electron, (right) U|| is the radiated spectral power polarized 
parallel to the magnetic field. 



which is similar to the parameters of a small synchrotron. The two polarized 
components of the radiation are u_l = fioe^c/lGn^u'^ (Jft^SxY + '^{Sx)'^) and 
U|, = /zoe^c/ievr^cu^ (s^^i^Sy) cos9 -^{S,) sm9f+{Q (Sy) cos ^ - 53 (S,) sin^)^ 
i.e. the components perpendicular and parallel to the plane of observation. 
Such calculations can be performed in a matter of minutes on a single pro- 
cessor (the figure shows 500 frequency bins by 500 angular bins resolution 
and took 30 minutes to calculate on a single processor of a 2 x 1.4 GHz iMac, 
but a calculation of 100 frequency bins by 100 angular bins resolution took 
~1 minute, including the calculation of the particle trajectory.) 

3.2. Linear and non-linear Thomson scattering 

To test non-relativistic motion extending to moderately relativistic mo- 
tion, Thomson scattering of a laser pulse from a single electron can be ana- 
lyzed. In non-relativistic Thomson scattering from a linearly polarized plane 
wave, an analytic solution for the radiated power is possible. The field 
strength corresponding to the interaction is parameterized by the normal- 
ized peak vector potential Oq = eEo/mciJo, where £^0 is the peak electric field 
strength and uq is the angular frequency of the laser. Returning to equation 
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Figure 4: Polar plot of radiation emitted by an electron oscillating in a linearly polarized 
plane electromagnetic wave, with low electric field strength (having a normalized vector 
potential of oq = 0.01) with the field vector along the x axis, calculated numerically. The 
vertical axis is the frequency normalized to the laser frequency, and the horizontal plane 
gives the radiated intensity. The angular shape of the structure is a set of sin^ Q lobes, in 
agreement with the analytic expression. 



|T|for a single particle, and inserting a time varying velocity (3 = ao sin(co'o^), 
and using a coordinate system so that angle B is measured between the 
polarization axis and the observation direction, the spectral power is: 



(fl /ioc^c 2 
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dudVt IGtt'^ 

To first-order in oq, assuming ag ^ 1, this can be integrated to give: 



(31) 
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Here, the Dirac delta function 6{x) is a representation of (A;/7r)sinc^(a;/fc) 
in the limit that k oo. The radiated energy is therefore distributed at a 
single frequency with an angular structure consisting of a sin^ shape. In 
figure |4| the angular distribution of spectral power is shown from calcula- 
tions using the quadratic algorithm with a time-step of cuqAt = 7r/25, with 
a plane electromagnetic wave polarized along the x-axis with a field strength 
of Oo = 0.01. Time is integrated to include 20 wave-periods. The angular 
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shape of the structure is a set of sin^ B lobes, in agreement with the ana- 
lytic expression. Although the spectrum as a function of energy is a peaked 
distribution centered at u = ujq, it is not a Dirac delta function. However, 
the analytic solution is for an infinite summation of wave-periods, whereas 
for obvious reasons a finite summation is calculated numerically. In terms 
of frequency, the calculation yields a sinc^(co' — ujq) shape, characteristic of a 
finite window, which approaches the Dirac delta as the size of the window is 
increased. 

As the field strength increases in the interaction, the electron motion 
becomes more complicated due to the magnetic field, and tends towards a 
figure of eight motion. This means that additional harmonics of the motion 
appear in the spectrum. For the case of exact back-scatter from an electron 
with an initial Lorentz factor of 70, the motion of the electron means that 
it experiences a Doppler shifted plane wave, and therefore the period of 
oscillation is modified. In reference [23] an analytic expression for the shifted 
fundamental frequency, uji, was given as: 

where (3o is the initial velocity of the particle, and the power per unit solid 
angle in the mth harmonic of uJi is: 



^ J 16. V7o^lSf ) S [^(™-l)/2(^C) - J(m+l)/2(mC)]' m odd , ^^^^ 

I m even , 

where Ju{x) is a Bessel function of the first kind and ( = al/{2al + 4). In 
figure [5} the power in harmonics of the fundamental frequency, normalized to 
the first harmonic, as a function of harmonic m, due to non-linear Thomson 
backscattering from a 5 MeV electron beam colliding with 10 periods of a 
linearly polarized plane wave with a normalized vector potential of oq = 
1 is shown. The temporal resolution is such that ujqAt = 7r/50, which 
represents only 100 sample points to calculate the spectrum from. Also shown 
is the (normalized) power from the analytic solution. Again the analytic 
solution models an infinite plane wave solution, and hence the harmonics are 
discrete. In the numerical solution, a finite number of periods is calculated, 
and therefore the harmonics have a sinc^(ct; — Um) shape. The scaling of the 
amplitude of the peaks and the position of the harmonics agrees well between 
the two solutions. For the numerical solution, the spectrum was calculated 
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Figure 5: The power in harmonics of the fundamental frequency, normahzed to the first 
harmonic, due to non-hnear Thomson backscattering from a 5 MeV electron beam colliding 
with a linear polarized plane wave with a normalized vector potential of ag = 1. (red) 
Calculated from the analytic expression of reference [21]. (blue) Calculated from numerical 
algorithm. 



as a function of frequency and then divided by the analytically calculated 
frequency Ui to give the horizontal axis. The 12th harmonic corresponds to 
a frequency of 784c(Jo, which for an 800 nm laser interaction corresponds to 
a photon energy of 1.2 keV. 

Figure |6] shows the angularly resolved plot of the same calculations as 
figure [sj and can be compared with figure 2(b) in reference [25], which shows 
a similar plot calculated from an analytic solution. Note that in their figure, 
only the first three harmonics were calculated. In the top right of the mesh- 
plot in figure [6} the next 3 harmonics can be observed. The finite width of the 
pulse, in both the numerical calculations here and the analytic calculations 
of reference [25|, results in the sinc^(co' — Um) shape of the spectral peaks. 

It should be noted that the radiation reaction force model of reference 
|26j has also been included to particle motions for other laser-particle inter- 
actions, and the energy loss by the particle compared with the energy in the 
calculated spectra, with good agreement. The question arises as to what reso- 
lution is required to accurately reproduce spectra using this method. Firstly, 
the particle trajectories themselves have to be reproduced accurately. Sec- 
ondly, the quadratic interpolation must be an accurate representation of the 
function KaX". Since the calculation is an integral, slight discontinuities at 
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Figure 6: The normalized spectral intensity as a function of normalized frequency, u jA'y^ujQ 
and angle 70^0 of radiation scattered by a 5 MeV electron beam colliding with a linear 
polarized plane wave with a normalized vector potential of ag = 1. (c.f. reference |25j ) 



grid points are not as important as the accuracy of the function between 
grid-points. This means that a sphne interpolation method is not necessar- 
ily better, and has been found to be worse, than an interpolation that only 
relates to local grid-points. If the coefficient for a cubic interpolation term 
is small, XsnT^ ^ 1, then the spectrum should be accurately modeled. The 
result of this is that the algorithm is not particularly efficient for calculations 
of highly relativistic particles (^GeV) performing large radius orbits such 
as in a classic synchrotron, although it will yield correct results given suffi- 
cient resolution (the required resolution for such orbits scales as 7, because 
the trigonometric nature of the motion means a cubic term exists and 
scales with 7). For such calculations other methods may be preferred. How- 
ever, for laser or plasma based accelerator interactions, for example, where 
typically the electron energies are < 10 GeV but the oscillation frequencies 
are very fast, 10^'^ — 10^^ Hz, the algorithm works well. 
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